//-------------------------------------------------------------------------------------------//
// Copyright(c) Wang Yong
// Author     : Wang Yong
// Date       : 2021-06-01
// Description: Refine the triangle mesh by split the marked element. It is designed for
// MoM adaptive mesh function.
//-------------------------------------------------------------------------------------------//
#include <cmath>
#include <algorithm>
#include "RefineTria.h"

// find all common edge
void RefineTria::fn_FindTriaEdges(std::vector<Node> &vR, std::vector<Tria> &vElem)
{
    constexpr int LocNode[3][2] = { {1,2}, {2,0}, {0,1} };

    if (vElem.empty()) return;

    int iNode, m, n, j, k;
    bool bFind;

    m_vTriaEdge.resize(vElem.size(), std::vector<int>(3));
    auto vvNode2Edge = std::vector<std::vector<int>>(vR.size());
    // set first tetrahedron
    std::vector<std::vector<int>> vvLib;
    auto vEdge = std::vector<int>(2);
    for (auto i = 0; i < 3; ++i) {
        m_vTriaEdge[0][i] = i;
        for (j = 0; j < 2; ++j) {
            iNode = vElem[0].node[LocNode[i][j]];
            vEdge[j] = iNode;
            vvNode2Edge[iNode].push_back(i);
        }
        std::sort(vEdge.begin(), vEdge.end());
        vvLib.push_back(vEdge);
    }
    // check the other tetrahedron
    for (auto i = 1; i < vElem.size(); ++i) {
        // check each edge
        for (n = 0; n < 3; ++n) {
            vEdge[0] = vElem[i].node[LocNode[n][0]];
            vEdge[1] = vElem[i].node[LocNode[n][1]];
            std::sort(vEdge.begin(), vEdge.end());
            bFind = false;
            for (j = 0; j < 2; ++j) {
                iNode = vEdge[j];
                for (k = 0; k < vvNode2Edge[iNode].size(); ++k) {
                    m = vvNode2Edge[iNode][k];
                    if (vvLib[m][0] == vEdge[0] && vvLib[m][1] == vEdge[1]) {
                        bFind = true;
                        m_vTriaEdge[i][n] = m;
                        break;
                    }
                }
                if (bFind) break;
            }
            if (!bFind) {
                m = int(vvLib.size());
                m_vTriaEdge[i][n] = m;
                vvNode2Edge[vEdge[0]].push_back(m);
                vvNode2Edge[vEdge[1]].push_back(m);
                vvLib.push_back(vEdge);
            }
        }
    }
    // get each edge length
    Node r;
    Real dist;
    m_vEdgeLen.reserve(vvLib.size());
    for (auto& v : vvLib) {
        r.x = vR[v[0]].x - vR[v[1]].x;
        r.y = vR[v[0]].y - vR[v[1]].y;
        r.z = vR[v[0]].z - vR[v[1]].z;
        dist = sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
        m_vEdgeLen.push_back(dist);
    }
}

// mark the edges need to be splitted
void RefineTria::fn_MarkSplitEdge(std::vector<Tria> &vElem, std::vector<bool> &vRefine)
{
    if (vElem.empty()) return;

    Real dMinEdge, dMaxEdge;
    int iMaxEdge, j;

    m_vEdgeFlag.resize(m_vEdgeLen.size(), false);
    // check all element edge
    for (auto i = 0; i < vElem.size(); ++i) {
        if (vRefine[i]) {
            iMaxEdge = 0;
            dMinEdge = dMaxEdge = m_vEdgeLen[m_vTriaEdge[i][0]];
            for (j = 1; j < 3; ++j) {
                if (dMinEdge > m_vEdgeLen[m_vTriaEdge[i][j]]) {
                    dMinEdge = m_vEdgeLen[m_vTriaEdge[i][j]];
                }
                if (dMaxEdge < m_vEdgeLen[m_vTriaEdge[i][j]]) {
                    dMaxEdge = m_vEdgeLen[m_vTriaEdge[i][j]];
                    iMaxEdge = j;
                }
            }
            dMinEdge *= 1.2;

            m_vEdgeFlag[m_vTriaEdge[i][iMaxEdge]] = true;
            for (j = 0; j < 3; ++j) {
                if (m_vEdgeLen[m_vTriaEdge[i][j]] > dMinEdge) {
                    m_vEdgeFlag[m_vTriaEdge[i][j]] = true;
                }
            }
        }
    }
}

// adapt the edge to improve the mesh quality
void RefineTria::fn_AdaptTriaEdge(std::vector<Tria> &vElem)
{
    Real dMinEdge;
    int edge[3], iEdge, m, j;

    for (auto i = 0; i < vElem.size(); ++i) {
        m = 0;
        for (j = 0; j < 3; ++j) {
            if (m_vEdgeFlag[m_vTriaEdge[i][j]]) {
                edge[m] = m_vTriaEdge[i][j];
                ++m;
            }
        }
        if (m > 0) {
            dMinEdge = m_vEdgeLen[edge[0]];
            for (j = 1; j < m; ++j) {
                if (dMinEdge > m_vEdgeLen[edge[j]]) {
                    dMinEdge = m_vEdgeLen[edge[j]];
                }
            }
            for (j = 0; j < 3; ++j) {
                iEdge = m_vTriaEdge[i][j];
                if (m_vEdgeLen[iEdge] > dMinEdge) {
                    m_vEdgeFlag[iEdge] = true;
                }
            }
        }
    }
}

// insert new point on each edge
void RefineTria::fn_AddEdgesPoint(std::vector<Tria> &vElem, std::vector<Node> &vR)
{
    constexpr int LocNode[3][2] = { {1,2}, {2,0}, {0,1} };

    if (vElem.empty()) return;

    Node r;
    int iEdge, j;

    // count how many point will be created;
    iEdge = 0;
    for (auto cut : m_vEdgeFlag) {
        if (cut) ++iEdge;
    }
    vR.reserve(vR.size() + iEdge);

    m_vEdgeNode.resize(m_vEdgeLen.size(), -1);
    // check all element edge
    for (auto i = 0; i < vElem.size(); ++i) {
        for (j = 0; j < 3; ++j) {
            iEdge = m_vTriaEdge[i][j];
            if (m_vEdgeFlag[iEdge] && m_vEdgeNode[iEdge] < 0) {
                m_vEdgeNode[iEdge] = int(vR.size());
                r.x = vR[vElem[i].node[LocNode[j][0]]].x + vR[vElem[i].node[LocNode[j][1]]].x;
                r.y = vR[vElem[i].node[LocNode[j][0]]].y + vR[vElem[i].node[LocNode[j][1]]].y;
                r.z = vR[vElem[i].node[LocNode[j][0]]].z + vR[vElem[i].node[LocNode[j][1]]].z;
                r.x *= 0.5;
                r.y *= 0.5;
                r.z *= 0.5;
                vR.push_back(r);
            }
        }
    }
}

// split each element
void RefineTria::fn_SplitEachTria(std::vector<Tria> &vElem)
{
    if (vElem.empty()) return;

    std::vector<Tria> vNewElem;

    int iCutEdge, j;
    for (auto i = 0; i < vElem.size(); ++i) {
        // count how many edges need to be splitted
        iCutEdge = 0;
        for (j = 0; j < 3; ++j) {
            if (m_vEdgeFlag[m_vTriaEdge[i][j]]) ++iCutEdge;
        }
 
        if (iCutEdge == 3) {
            // split all three triangle edges
            fn_SplitEdge3(vElem[i], m_vTriaEdge[i], vNewElem);
        } else if (iCutEdge == 2) {
            // split two triangle edges
            fn_SplitEdge2(vElem[i], m_vTriaEdge[i], vNewElem);
        } else if (iCutEdge == 1) {
            // split only one triangle edges
            fn_SplitEdge1(vElem[i], m_vTriaEdge[i], vNewElem);
        } else {
            vNewElem.push_back(vElem[i]);
        }
    }

    vElem.clear();
    vNewElem.shrink_to_fit();
    vElem = vNewElem;
}

// split only one triangle edges
void RefineTria::fn_SplitEdge1(Tria &oldTria, std::vector<int> &vEdgeIndx, std::vector<Tria> &vNewTria)
{
    Tria tri;

    tri.id = oldTria.id;

    if (m_vEdgeFlag[vEdgeIndx[0]]) {
        // create the first triangle
        tri.node[0] = oldTria.node[0];
        tri.node[1] = oldTria.node[1];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
        vNewTria.push_back(tri);

        // create the second triangle
        tri.node[0] = oldTria.node[0];
        tri.node[1] = m_vEdgeNode[vEdgeIndx[0]];
        tri.node[2] = oldTria.node[2];
        vNewTria.push_back(tri);
    } else if (m_vEdgeFlag[vEdgeIndx[1]]) {
        // create the first triangle
        tri.node[0] = oldTria.node[0];
        tri.node[1] = oldTria.node[1];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
        vNewTria.push_back(tri);

        // create the second triangle
        tri.node[0] = oldTria.node[1];
        tri.node[1] = oldTria.node[2];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
        vNewTria.push_back(tri);
    } else {
        // create the first triangle
        tri.node[0] = oldTria.node[0];
        tri.node[1] = m_vEdgeNode[vEdgeIndx[2]];
        tri.node[2] = oldTria.node[2];
        vNewTria.push_back(tri);

        // create the second triangle
        tri.node[0] = oldTria.node[1];
        tri.node[1] = oldTria.node[2];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
        vNewTria.push_back(tri);
    }
}

// split two triangle edges
void RefineTria::fn_SplitEdge2(Tria &oldTria, std::vector<int> &vEdgeIndx, std::vector<Tria> &vNewTria)
{
    Tria tri;

    tri.id = oldTria.id;

    if (!m_vEdgeFlag[vEdgeIndx[0]]) {
        // create the first triangle
        tri.node[0] = oldTria.node[0];
        tri.node[1] = m_vEdgeNode[vEdgeIndx[2]];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
        vNewTria.push_back(tri);

        // create the second and third triangle
        if (m_vEdgeLen[vEdgeIndx[1]] > m_vEdgeLen[vEdgeIndx[2]]) {
            tri.node[0] = oldTria.node[1];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[1]];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
            vNewTria.push_back(tri);

            tri.node[0] = oldTria.node[1];
            tri.node[1] = oldTria.node[2];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
            vNewTria.push_back(tri);
        } else {
            tri.node[0] = oldTria.node[2];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[1]];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
            vNewTria.push_back(tri);

            tri.node[0] = oldTria.node[1];
            tri.node[1] = oldTria.node[2];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
            vNewTria.push_back(tri);
        }
    } else if (!m_vEdgeFlag[vEdgeIndx[1]]) {
        // create the first triangle
        tri.node[0] = oldTria.node[1];
        tri.node[1] = m_vEdgeNode[vEdgeIndx[0]];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
        vNewTria.push_back(tri);

        // create the second and third triangle
        if (m_vEdgeLen[vEdgeIndx[0]] > m_vEdgeLen[vEdgeIndx[2]]) {
            tri.node[0] = oldTria.node[0];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[2]];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
            vNewTria.push_back(tri);

            tri.node[0] = oldTria.node[2];
            tri.node[1] = oldTria.node[0];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
            vNewTria.push_back(tri);
        } else {
            tri.node[0] = oldTria.node[0];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[2]];
            tri.node[2] = oldTria.node[2];
            vNewTria.push_back(tri);

            tri.node[0] = m_vEdgeNode[vEdgeIndx[0]];
            tri.node[1] = oldTria.node[2];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
            vNewTria.push_back(tri);
        }
    } else {
        // create the first triangle
        tri.node[0] = oldTria.node[2];
        tri.node[1] = m_vEdgeNode[vEdgeIndx[1]];
        tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
        vNewTria.push_back(tri);

        // create the second and third triangle
        if (m_vEdgeLen[vEdgeIndx[0]] > m_vEdgeLen[vEdgeIndx[1]]) {
            tri.node[0] = oldTria.node[0];
            tri.node[1] = oldTria.node[1];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
            vNewTria.push_back(tri);

            tri.node[0] = oldTria.node[0];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[0]];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
            vNewTria.push_back(tri);
        } else {
            tri.node[0] = oldTria.node[1];
            tri.node[1] = m_vEdgeNode[vEdgeIndx[0]];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
            vNewTria.push_back(tri);

            tri.node[0] = oldTria.node[0];
            tri.node[1] = oldTria.node[1];
            tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
            vNewTria.push_back(tri);
        }
    }
}

// split all three triangle edges
void RefineTria::fn_SplitEdge3(Tria &oldTria, std::vector<int> &vEdgeIndx, std::vector<Tria> &vNewTria)
{
    Tria tri;

    tri.id = oldTria.id;
    // create the first triangle
    tri.node[0] = oldTria.node[0];
    tri.node[1] = m_vEdgeNode[vEdgeIndx[2]];
    tri.node[2] = m_vEdgeNode[vEdgeIndx[1]];
    vNewTria.push_back(tri);

    // create the second triangle
    tri.node[0] = oldTria.node[1];
    tri.node[1] = m_vEdgeNode[vEdgeIndx[0]];
    tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
    vNewTria.push_back(tri);

    // create the third triangle
    tri.node[0] = m_vEdgeNode[vEdgeIndx[0]];
    tri.node[1] = m_vEdgeNode[vEdgeIndx[1]];
    tri.node[2] = m_vEdgeNode[vEdgeIndx[2]];
    vNewTria.push_back(tri);

    // create the fourth triangle
    tri.node[0] = oldTria.node[2];
    tri.node[1] = m_vEdgeNode[vEdgeIndx[1]];
    tri.node[2] = m_vEdgeNode[vEdgeIndx[0]];
    vNewTria.push_back(tri);
}